Replication of general RNAseq analysis of the response to Nivolumab immunotherapy in cases of advanced melanoma

1 Abstract

Tumor progression during immunotherapy treatment is a complex process that has not yet been fully characterized due to multiple dynamic factors that can influence it. For this reason, we decided to use data from BioProject PRJNA356761, which corresponds to patients with advanced melanoma who had experienced progression with ipilimumab or who had not received prior treatment with ipilimumab, both before and after initiating treatment with nivolumab. The tumors were examined using whole exome sequencing, transcriptome and/or T-cell receptor (TCR) analysis. The results obtained provide us with evidence that Nivolumab treatment has better efficacy in recognizing and targeting cancer cells, in our enrichment analysis we observed that the most prominent genes are associated with acute inflammatory response, acute phase response, humoral immune response, metabolic and catabolic processes. Which makes sense considering the mechanisms of function of immunotherapies. In the bioproject article they obtained similar results to ours, they reported highly expressed genes related to the immune system, which makes sense considering the mechanisms of function of immunotherapies. This suggests that the treatment is effective and confirms the replicability of the analysis. The data provide information on the mechanism of action of nivolumab, however we consider that further studies could be performed and the sample size could be increased for a better understanding of the interactions of this treatment.

2 Introduction

Melanoma is a type of skin cancer that originates in pigment-producing cells called melanocytes. Immunotherapy has been positioned as one of the most effective treatments for this cancer, however access to these treatments is very limited due to cost. In addition, some genomic profiles have been observed that may favor or affect the response to treatments. In this case, we are analyzing data from patients who were treated with immune checkpoints and for this type of immunotherapy the mechanisms that modulate tumor evolution are not clear.

The genomic characteristics of a tumor may contribute to the response to checkpoint blockade (Li, et al., 2022). High tumor mutation loads increase the likelihood of generating immunogenic neoantigens, which facilitates the recognition of a tumor as foreign and thus having a higher number of clonal neoantigens elicits more effective immune responses (McGranahan, et al., 2016). Characteristics of the tumor microenvironment have also been linked to the response to checkpoint inhibitor therapy, primarily CD8+ T-cell infiltration.

Understanding how activation of the immune system, through checkpoint inhibition, affects tumor mutation and the tumor microenvironment (TME). In order to analyze genomic changes that could contribute to the understanding of these dynamics, we conducted a study of melanoma samples before and after treatment with nivolumab (Nivo; an anti-PD-1 agent).

3 Methods

First of all, a new directory in the DNA High-Performance Computing Cluster was created, which is

/mnt/Timina/bioinfoII/amarin/RNAseq2/

This directory has the next organization:

|-data/           # Directory to download the raw data
|-data_trimmed/   # Directory to save the post-QC data
|-kallisto_out/   # Directory to save the kallisto outputs
|-QC/             # Directory to perform raw QC
|-QC_trimmed/     # Directory to perform QC on the trimmed data
|-scripts/        # Some scripts used to send jobs

3.1 Download

The dataset choose to analysis has de GEO Accession GSE91061. Since downloading from the NCBI may be difficult and quite technical, it was choosen to use the European mirror, the European Nucleotide Archive from the EMBL-EBI. By browing on the GEO web page, we see that the data are from the BioProject PRJNA356761. We use that ID to search on ENA.

Due the fact that the are a lot of file needed to be downloaded, we use the web application SRA Explorer to get a script to get the bulk download.

We then select all the files and click on “Add to collection”. After that, we see our saves datasets and select the bash script for downloading FastQ files option. We save the script with the name downloadScript.sh, and save it in the HPCC, in the //data/ directory.

Then, we send a job using the script. It lasted around 5 days.

3.2 Raw Quality Control

Once the data is ready, we perform quality control. We just send a job using fastqc with all the fastq files via a wildcard:

fastqc /mnt/Timina/bioinfoII/amarin/RNAseq2/data/*fastq.gz -o /mnt/Timina/bioinfoII/amarin/RNAseq2/QC

Then, we used multiqc to visualize all the samples at once:

multiqc /mnt/Timina/bioinfoII/amarin/RNAseq2/QC -o /mnt/Timina/bioinfoII/amarin/RNAseq2/QC

This allowed us to see some irregularities. First, the first to nucleotides of each read had a smaller mean quality score, compared to the rest of the read.

Besides, in the heatmap, the base sequence sequence content seemed to be biased on the first base pairs of the read.

All the samples resembled more or less the next percentage of base content.

Nevertheless, MultiQC didn’t detect any adapter content. Based on that, we decided to trimm the first 12 nucleotides of each pair.

3.3 Trimming

We send a job in the //data/ directory with a script using the software Trimmomatic:

for i in *_1.fastq.gz;
    do echo
    trimmomatic PE -threads 40 -phred33 -trimlog trim.log \
        $i "${i%_1.fastq.gz}_2.fastq.gz" \
        ../data_trimmed/"${i%_1.fastq.gz}_1_trimmed.fastq.gz" \
        ../data_trimmed/"${i%_1.fastq.gz}_1_unpaired.fastq.gz" \
        ../data_trimmed/"${i%_1.fastq.gz}_2_trimmed.fastq.gz" \
        ../data_trimmed/"${i%_1.fastq.gz}_2_unpaired.fastq.gz" \
        HEADCROP:12
  done

Here we just cut the first 12 nucleotides of each read and create a trim.log log file by the use of 40 threads. We know that the Phred score are based on ASCII 33 because previously we found a # symbol on the reads.

This job took a couple of days.

3.4 Trimmed Quality Control

Once the data was trimmed, we did one more time the quality control.

fastqc /mnt/Timina/bioinfoII/amarin/RNAseq2/data_trimmed/*fastq.gz \
  -o /mnt/Timina/bioinfoII/amarin/RNAseq2/QC_trimmed
multiqc /mnt/Timina/bioinfoII/amarin/RNAseq2/QC_trimmed \
  -o /mnt/Timina/bioinfoII/amarin/RNAseq2/QC_trimmed

This time, the quality of the first pair bases was better.

Besides, the “per base sequence content” was less biased:

Each sample was more homogeneous.

Besides, other indicators, such as the “per base N content” and “per sequence GC content” performed better.

3.5 Pseudoaligment to reference genome

Once the data was ready, we decided to generate count matrices using Kallisto. To perform this, we first need a reference transcriptome. We choose the latest version available at GENCODE, which is the release 43. The file was downloaded at the parent directory of this project:

wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_43/gencode.v43.transcripts.fa.gz

Then, generated the kallisto index:

kallisto index -i /mnt/Timina/bioinfoII/amarin/RNAseq2/kallisto_out/Hs_ref.kidx \
  /mnt/Timina/bioinfoII/amarin/RNAseq2/gencode.v43.transcripts.fa.gz

Then, we used a bash script to perform the quantification with all the samples.

for file in /mnt/Timina/bioinfoII/amarin/RNAseq2/data_trimmed/*_1_trimmed.fastq.gz
do
    op=$(basename $file | cut -c-10 )         
    file2=$(echo $file | sed 's/_1_/_2_/')  
    mkdir /mnt/Timina/bioinfoII/amarin/RNAseq2/kallisto_out/${op}
    
    kallisto quant -i /mnt/Timina/bioinfoII/amarin/RNAseq2/kallisto_out/Hs_ref.kidx \
    -o /mnt/Timina/bioinfoII/amarin/RNAseq2/kallisto_out/${op} \
    --threads 40 $file $file2
done

3.6 Differencial Expresion

For the differential expression analysis we used DESeq2, a Bioconductor package running on R 4.3.0 using a PC. The next, are the libraries to use.

library(plotly)
library(tximport)
library(tidyverse)
library(DESeq2)
library(ggplot2)
library(ggrepel) 
library(rhdf5)
library(GenomicFeatures)
library(dplyr)
library(biomaRt)
library(gprofiler2)
library(ComplexHeatmap)
library(genefilter)
library(clusterProfiler)
library(AnnotationDbi)
library(org.Hs.eg.db)
#library(ggupset)
library(enrichplot)

3.6.1 Data importation

First, we imported the kallisto quantification tables to a local computer from the DNA HPCC. Then, we created a tx2gene table, which maps annotated genes with transcripts. To do this, we used the Bioconductor package biomaRt.

ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")


tx2gene = getBM(attributes=c('ensembl_transcript_id','ensembl_gene_id'), 
                   mart = ensembl)

Next, we import the quantification table to the R environment. To do this, we use the Bioconductor package rhdf5.

## Paths to h5 archives
files = file.path(list.dirs("./kallisto_data", recursive = FALSE),
                  "abundance.h5")
names(files) = str_extract(files, "SRR\\d+")

## Importation
txi.kallisto.tsv = tximport(files, type = "kallisto", tx2gene = tx2gene,
                            ignoreAfterBar = TRUE, ignoreTxVersion=TRUE)

## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 
## summarizing abundance
## summarizing counts
## summarizing length

Next, we create a metadata table to know the treatment for each SRA sample. Since it wasn’t available on SRA Run Selector web page, we needed to manually create it.

## Metadata importation to R
meta = read.csv("./metadata.csv", header = FALSE)
samples = column_to_rownames(meta, var = "V1")

Now, we create a DESeq2 dataset-like object, called DESeqDataSet, using the imported counts, sample names and metadata table.

dds = DESeqDataSetFromTximport(txi.kallisto.tsv,
                               colData = samples,
                               design = ~V2)

## using counts and average transcript lengths from tximport

3.6.2 Exploratory Analysis

Before doing some exploration analysis, we pre-filter low counts.

## Pre-filtering

keep = rowSums(counts(dds)) >= 10
dds = dds[keep, ]

## Factor level
dds$V2 = factor(dds$V2, levels = c("Pre", "On"))

To identify if there’s a batch effect, we graph a PCA; normalizing the data with the variance stabilizing transformation.

vst = vst(dds, blind = FALSE)

Once we hve the normalized data, we graph a tridimensional and interactive PCA plot.

pca = prcomp(t(assay(vst)), rank. = 3)

components = pca[["x"]]
components = data.frame(components)

var = summary(pca)[["importance"]]['Proportion of Variance',]
var = 100*sum(var)

fig = plot_ly(components, x = ~PC1, y = ~PC2, z = ~PC3, color = dds$V2)
fig = fig %>% layout(title = "PCA")
fig

As we can observe, there is not a visible cluster suggesting a batch effect. In fact, by some reason, almos all the samples cluster into a big cloud, regardless to their treatment.

3.6.3 Differential Expression Analysis

Once we know there are not batch effects in our data, we perform the Differential Expression Analysis.

dds = DESeq(dds)

To consider significant genes, we will use a significance level of \(\alpha=0.01\), and a log fold change expression greater than 1 or lower than -1.

res = results(dds, alpha = 0.01)

## Significative genes
resSig = subset(res, padj < 0.01)
dim(resSig)

## [1] 80  6

We have 80 significant genes. We can visualize them in a volcano plot.

res.df = as.data.frame(res)

## Add genes symbols
res.df$symbols = mapIds(org.Hs.eg.db, keys = rownames(res.df),
                        keytype = "ENSEMBL", column = "SYMBOL")


### Up and down genes
res.df$diffexpressed = "NO"
res.df$diffexpressed[res.df$log2FoldChange > 1 & res.df$padj < 0.01] <- "UP"
res.df$diffexpressed[res.df$log2FoldChange < -1 & res.df$padj < 0.01] <- "DOWN"



ggplot(data = res.df, aes(x=log2FoldChange, y=-log10(padj), col=diffexpressed, label=symbols)) + 
  geom_point() 

3.6.4 Significant genes extraction

Once we have identified the significant genes, we extract them in tsv files. One for the up-regulated ones, and another for the down-regulated ones.

## Up genes
resSigUP = subset(resSig, log2FoldChange >= 1)
resSigUP = cbind(rownames(resSigUP), resSigUP) # Converts rownames into first column
colnames(resSigUP)[1] = "GENE_ID" # Adds name to first column

write.table(resSigUP, file = "./DEG_results/UP_genes.tsv", quote = FALSE,
            sep = "\t",row.names = FALSE)


## Down genes
resSigDOWN = subset(resSig, log2FoldChange <= -1)
resSigDOWN = cbind(rownames(resSigDOWN), resSigDOWN)
colnames(resSigDOWN)[1] = "GENE_ID"

write.table(resSigDOWN, file = "./DEG_results/DOWN_genes.tsv", quote = FALSE,
            sep = "\t",row.names = FALSE)

We got 6 significant downregulated genes and 73 upregulated significant genes.

3.7 Gene Ontology Enrichment

To make the GO terms enrichment, we use the org.Hs.eg.db database, containing the annotated genes from Homo sapiens. The ontology was made using the Biological Process (BP).

genes = rownames(resSig[resSig$log2FoldChange >= 1, ])

GO = enrichGO(gene = genes,
              OrgDb = org.Hs.eg.db,
              keyType = "ENSEMBL", 
              ont = "BP")

GO.df = as.data.frame(GO)

3.7.1 Gene-Concept Network

We can visualize the enriched genes depending on their biological process.

GO2 = setReadable(GO, OrgDb = org.Hs.eg.db)

cnetplot(GO2, node_label = "all",color_category='firebrick', 
          color_gene='steelblue', showCategory = 10)

4 Results and Discussion

Our results provide evidence that the Immunotherapy with Nivolumab helps the immune system recognize and attack cancer cells more effectively. In our differential gene expression analysis we found that 6 of our significant genes were downregulated and 73 were upregulated in patients being treated with Nivolumab. We then performed a Gene Ontology (GO) enrichment analysis to gain a better understanding of the functional implications of our genes. This analysis show us that the enriched genes that we found were associated with the acute inflammatory response, acute-phase response, humoral immune response, metabolic and catabolic processes which makes sense as immunotherapies work by targeting specific proteins on immune cells and cancer cells that regulate immune responses.

Cancerous tumor cells have the ability to produce signals that cause inflammation. Cancer can be affected by inflammation in both favorable and unfavorable ways. On one hand, it might help to attract immune cells like neutrophils and macrophages that can recognize and eliminate cancer cells.On the other hand, chronic inflammation can create a tumor-promoting environment by promoting angiogenesis and tissue remodeling, which can facilitate tumor growth and metastasis (Li, L. et al., 2020).

The tumor or the host immune system’s response to the tumor can both initiate the acute-phase response in cancer (Hugo Becerril Gonzalez et al., 2018). Acute-phase protein levels have been found to be higher in a variety of cancers. Acute-phase proteins can be monitored as diagnostic or prognostic indicators in some malignancies.

The humoral immune response can be modulated by tumor cells to evade immune surveillance. Tumors may downregulate the expression of antigens recognized by antibodies or produce factors that inhibit antibody production (Anushruti Sarvaria et al., 2017). However, in some cases, the humoral immune response can contribute to anti-tumor immunity and be harnessed for therapeutic purposes.

In the the paper of the bioproject they mentioned that highly expressed genes were immune-related, suggesting pre-existing immune recognition of the tumor. By comparing our results with the paper we noticed that they had similar results as they also found differentially expressed genes related to host immunity, which increased their expression and with metabolism, which decreased their expression in responders. With this comparison, we can be sure that our analysis was perform correctly and that our results are assuring the efficacy of the immunotherapy.

5 Conclusions

In summary, the study showed that immunotherapy with Nivolumab improves the immune system’s ability to target and fight cancer cells. Gene expression analysis revealed changes in immune-related genes, inflammation, and metabolic processes. The findings support the effectiveness of Nivolumab and provide insights into its mechanism of action. Overall, the importance of these findings lies in advancing our knowledge of immunotherapy, providing insights into the molecular mechanisms underlying its effectiveness, and potentially guiding future research and clinical approaches in cancer treatment.

6 References

Li, L., Yu, R., Cai, T., Chen, Z., Lan, M., Zou, T., Wang, B., Wang, Q., Zhao, Y., & Cai, Y. (2020). Effects of immune cells and cytokines on inflammation and immunosuppression in the tumor microenvironment. 88, 106939–106939. https://doi.org/10.1016/j.intimp.2020.106939

Zhao, H., Wu, L., Yan, G., Chen, Y., Zhou, M., Wu, Y., & Li, Y. (2021). Inflammation and tumor progression: signaling pathways and targeted intervention. 6(1). https://doi.org/10.1038/s41392-021-00658-5

Hugo Becerril Gonzalez, Catharina Hagerling, & Werb, Z. (2018). Roles of the immune system in cancer: from tumor initiation to metastatic progression. 32(19-20), 1267–1284. https://doi.org/10.1101/gad.314617.118

Anushruti Sarvaria, J. Alejandro Madrigal, & Saudemont, A. (2017). B cell regulation in cancer and anti-tumor immunity. 14(8), 662–674. https://doi.org/10.1038/cmi.2017.35

Li, Y., Wang, Q., Chen, Y. et al. Somatic mutations in DCC are associated with genomic instability and favourable outcomes in melanoma patients treated with immune checkpoint inhibitors. Br J Cancer 127, 1411–1423 (2022). https://doi.org/10.1038/s41416-022-01921-4

McGranahan N, Furness AJ, Rosenthal R, Ramskov S, Lyngaa R, Saini SK, Jamal-Hanjani M, Wilson GA, Birkbak NJ, Hiley CT, Watkins TB, Shafi S, Murugaesu N, Mitter R, Akarca AU, Linares J, Marafioti T, Henry JY, Van Allen EM, Miao D, Schilling B, Schadendorf D, Garraway LA, Makarov V, Rizvi NA, Snyder A, Hellmann MD, Merghoub T, Wolchok JD, Shukla SA, Wu CJ, Peggs KS, Chan TA, Hadrup SR, Quezada SA, Swanton C. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science. 2016 Mar 25;351(6280):1463-9. doi: 10.1126/science.aaf1490. Epub 2016 Mar 3. PMID: 26940869; PMCID: PMC4984254.

7 Session Information

sessionInfo()

## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Pop!_OS 22.04 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/Mexico_City
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] enrichplot_1.20.0           org.Hs.eg.db_3.17.0        
##  [3] clusterProfiler_4.8.1       genefilter_1.82.1          
##  [5] ComplexHeatmap_2.16.0       gprofiler2_0.2.1           
##  [7] biomaRt_2.56.0              GenomicFeatures_1.52.0     
##  [9] AnnotationDbi_1.62.1        rhdf5_2.44.0               
## [11] ggrepel_0.9.3               DESeq2_1.40.1              
## [13] SummarizedExperiment_1.30.1 Biobase_2.60.0             
## [15] MatrixGenerics_1.12.0       matrixStats_0.63.0         
## [17] GenomicRanges_1.52.0        GenomeInfoDb_1.36.0        
## [19] IRanges_2.34.0              S4Vectors_0.38.1           
## [21] BiocGenerics_0.46.0         lubridate_1.9.2            
## [23] forcats_1.0.0               stringr_1.5.0              
## [25] dplyr_1.1.0                 purrr_1.0.1                
## [27] readr_2.1.4                 tidyr_1.3.0                
## [29] tibble_3.1.8                tidyverse_2.0.0            
## [31] tximport_1.28.0             plotly_4.10.1              
## [33] ggplot2_3.4.1              
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.3.0            BiocIO_1.10.0            bitops_1.0-7            
##   [4] ggplotify_0.1.0          filelock_1.0.2           polyclip_1.10-4         
##   [7] XML_3.99-0.14            lifecycle_1.0.3          doParallel_1.0.17       
##  [10] lattice_0.21-8           MASS_7.3-59              crosstalk_1.2.0         
##  [13] magrittr_2.0.3           sass_0.4.5               rmarkdown_2.20          
##  [16] jquerylib_0.1.4          yaml_2.3.7               cowplot_1.1.1           
##  [19] DBI_1.1.3                RColorBrewer_1.1-3       zlibbioc_1.46.0         
##  [22] ggraph_2.1.0             RCurl_1.98-1.12          yulab.utils_0.0.6       
##  [25] tweenr_2.0.2             rappdirs_0.3.3           circlize_0.4.15         
##  [28] GenomeInfoDbData_1.2.10  tidytree_0.4.2           annotate_1.78.0         
##  [31] codetools_0.2-19         DelayedArray_0.26.2      DOSE_3.26.1             
##  [34] xml2_1.3.3               ggforce_0.4.1            tidyselect_1.2.0        
##  [37] shape_1.4.6              aplot_0.1.10             farver_2.1.1            
##  [40] viridis_0.6.3            BiocFileCache_2.8.0      GenomicAlignments_1.36.0
##  [43] jsonlite_1.8.4           GetoptLong_1.0.5         ellipsis_0.3.2          
##  [46] tidygraph_1.2.3          survival_3.5-5           iterators_1.0.14        
##  [49] foreach_1.5.2            tools_4.3.0              progress_1.2.2          
##  [52] treeio_1.24.0            Rcpp_1.0.10              glue_1.6.2              
##  [55] gridExtra_2.3            xfun_0.39                qvalue_2.32.0           
##  [58] withr_2.5.0              fastmap_1.1.0            rhdf5filters_1.12.1     
##  [61] fansi_1.0.4              digest_0.6.31            timechange_0.2.0        
##  [64] R6_2.5.1                 gridGraphics_0.5-1       colorspace_2.1-0        
##  [67] GO.db_3.17.0             RSQLite_2.3.1            utf8_1.2.3              
##  [70] generics_0.1.3           data.table_1.14.8        rtracklayer_1.60.0      
##  [73] prettyunits_1.1.1        graphlayouts_1.0.0       httr_1.4.4              
##  [76] htmlwidgets_1.6.2        S4Arrays_1.0.4           scatterpie_0.1.9        
##  [79] pkgconfig_2.0.3          gtable_0.3.1             blob_1.2.3              
##  [82] XVector_0.40.0           shadowtext_0.1.2         htmltools_0.5.4         
##  [85] fgsea_1.26.0             clue_0.3-64              scales_1.2.1            
##  [88] png_0.1-8                ggfun_0.0.9              knitr_1.42              
##  [91] rstudioapi_0.14          tzdb_0.3.0               reshape2_1.4.4          
##  [94] rjson_0.2.21             nlme_3.1-162             curl_5.0.0              
##  [97] cachem_1.0.6             GlobalOptions_0.1.2      parallel_4.3.0          
## [100] HDO.db_0.99.1            restfulr_0.0.15          pillar_1.8.1            
## [103] vctrs_0.5.2              dbplyr_2.3.0             xtable_1.8-4            
## [106] cluster_2.1.4            evaluate_0.20            cli_3.6.0               
## [109] locfit_1.5-9.7           compiler_4.3.0           Rsamtools_2.16.0        
## [112] rlang_1.0.6              crayon_1.5.2             labeling_0.4.2          
## [115] plyr_1.8.8               stringi_1.7.12           viridisLite_0.4.1       
## [118] BiocParallel_1.34.1      assertthat_0.2.1         munsell_0.5.0           
## [121] Biostrings_2.68.1        lazyeval_0.2.2           GOSemSim_2.26.0         
## [124] Matrix_1.5-1             patchwork_1.1.2          hms_1.1.2               
## [127] bit64_4.0.5              Rhdf5lib_1.22.0          KEGGREST_1.40.0         
## [130] highr_0.10               igraph_1.4.3             memoise_2.0.1           
## [133] bslib_0.4.2              ggtree_3.8.0             fastmatch_1.1-3         
## [136] bit_4.0.5                downloader_0.4           gson_0.1.0              
## [139] ape_5.7-1